home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / sixsig.exe / SSS.C next >
Text File  |  1992-07-11  |  33KB  |  997 lines
  1. /*
  2.  * Six Sigma Simulation (SSS) - This program graphically demonstrates
  3.  * histograms of pseudo-random i.i.d. Normal variables with mean zero and
  4.  * variance one.
  5.  *
  6.  * Source code for 5 functions are NOT included here because they are
  7.  * copyrighted materials that you may obtain directly from the publisher:
  8.  * ran1(), ran2(), ran3() and gasdev() are "portable" psuedo-random number
  9.  * generator functions of Press, et.al. (1988) "Numerical Recipes in C,"
  10.  * Ch.7, Cambridge University Press.  Similarly, erfcc() is the compact
  11.  * error function of Chapter 6 of that same book.
  12.  */
  13.  
  14. #include <graph.h>          // _outtext, _settextcolor, _settextposition
  15. #include <string.h>         // strlen
  16. #include <signal.h>         // signal handler
  17. #include <stdio.h>          // sprintf
  18. #include <stdlib.h>         // toupper
  19. #include <malloc.h>         // malloc, free
  20. #include <time.h>           // time, clock
  21. #include <conio.h>          // kbhit
  22. #include <math.h>
  23.  
  24. /* Macros */
  25. #define square( x )                   ( (x) * (x) )
  26. #define _outxcr( ach, col, row )      { _settextposition( row, col ); \
  27.                                         _outtext( ach ); }
  28. /* Constants */
  29. #define ESC        27
  30. #define BLANK      32
  31. #define BLOCK      223  /* top-half bar; full=219 & bottom-half=220 */
  32. #define EFRAC      0.75
  33. #define TOP        1
  34. #define LEFTCOLUMN 49
  35. #define PROMPTPOS  25
  36. #define WIDTH      31
  37. #define HEIGHT     (cszMenu + 2)
  38. #define SINGLCOLOR 7
  39. #define BACKCOLOR  0L
  40. #define MENUCOLOR  3
  41. #define MENUERASE  0
  42.  
  43. enum BOOL { FALSE, TRUE };
  44.  
  45. /* Structure type for colored bars */
  46.  
  47. typedef struct _HIST {
  48.     int   length;
  49.     short color;
  50.     } HIST;
  51.  
  52. /* Global variables */
  53.  
  54. char *aszMenu[] = {
  55.     "",
  56.     "     Six Sigma Simulator   ",
  57.     "",
  58.     "True",
  59.     "Cumulative",
  60.     "Incremental",
  61.     "Lines",
  62.     "MenuOff         NarrowCells",
  63.     "Parameters      ResetAccum ",
  64.     "SampleSize      WideCells  ",
  65.     "",
  66.     "",
  67.     "",
  68.     "",
  69.     "",
  70.     "Press a Menu Letter Key",
  71.     " ( T C I L M N P R S W ) ",
  72.     "or Press ESCape to quit: "
  73.     };
  74.  
  75. char   buf[256], yesno;
  76. char   savnam[13], parnam[13], Hbar[77], MenuCh[2], String[31];
  77. int    cszMenu = sizeof( aszMenu ) / sizeof( aszMenu[0] );
  78. int    dispMenu, dispType, topRow, botRow;
  79. int    reppbat, iter, tRows, iColorMax, dofinc, doftot, bwidth;
  80. int    idum, iduminit, firstobs, lastobs, firstcum, lastcum, prevC;
  81. int    bmtype, rntype, cellwid;
  82. time_t *ltime, timer;
  83. long   counts[50], totals[50], reps, maxreps;
  84. double zvalue, content, rroot2, pgreater, pequal;
  85. double chisqinc, chisqtot, exptmin, fwidth;
  86. double maxfrac, expfrac[50], exptail[50], cutval[50];
  87. double obsfrac[50], chisqobs[50], cumfrac[50], chisqcum[50];
  88. FILE   *savfile, *parfile;
  89. HIST   abar[50];              /* max histogram bars = rows on screen */
  90. struct videoconfig vc;
  91. struct tm *curtime;
  92.  
  93. /* Function prototypes */
  94. void   main( void );
  95. void   handler( int sig );
  96. void   ClrScreen( void );
  97. void   DrawCell( int iRow );
  98. void   DrawFrame( int iTop, int Left, int Width, int Height );
  99. void   DrawHist( int Code );
  100. void   DrawMenu( int Color );
  101. void   ObsrvHist( void );
  102. void   ResetHist( void );
  103. void   SSSparm( void );
  104. void   SampSize( void );
  105. int    RunMenu( void );
  106. double box_mull( int *idum );
  107. float  whran( int *idum );
  108. float  ranX( int *idum );
  109. int    initseed( void );
  110. time_t exptime( time_t *ltime );
  111.  
  112. extern void exit( int code );
  113. extern double erfcc( double x );
  114. extern double gasdev( int *idum );
  115. extern float  ran1( int *idum );
  116. extern float  ran2( int *idum );
  117. extern float  ran3( int *idum );
  118.  
  119. void handler( int sig ) {              /* Interrupt Handler */
  120.         signal( SIGINT, SIG_IGN );     /* ignore further signals */
  121.         _setvideomode( _DEFAULTMODE );
  122.         exit( 1 );
  123.         }
  124.  
  125. /*
  126.  * Brute Force Box-Muller Transformation.
  127.  */
  128. double box_mull( idum )
  129. int *idum; {
  130.  
  131.         static int iset=0;
  132.         static double next;
  133.         double rad,ang;
  134.  
  135.         if( iset == 0 ) {
  136.             rad = sqrt( -2.0 * log( (double)ranX( idum ) ) );
  137.             ang = 2.0 * 3.1415926536 * (double)ranX( idum );
  138.             next = rad * cos( ang );
  139.             iset = 1;
  140.             return( rad * sin( ang ) );
  141.             }
  142.         else {
  143.             iset=0;
  144.             return next;
  145.             }
  146.         }
  147.  
  148. void DrawMenu( int Color ) {    /* Color = 0 ==> erase */
  149.     int  i;
  150.     _settextcolor( Color );
  151.     _setbkcolor( BACKCOLOR );
  152.     DrawFrame( TOP, LEFTCOLUMN - 3, WIDTH + 3, HEIGHT );
  153.     for( i = 0; i < cszMenu; i++ )
  154.         _outxcr( aszMenu[i], LEFTCOLUMN, TOP + i + 1 );
  155.     }
  156.  
  157. /*
  158.  * DrawFrame - Makes a rectangular frame using the double-line box
  159.  * characters. The parameters iTop, iLeft, iWidth, and iHeight are
  160.  * the row and column arguments for the upper-left and lower-right
  161.  * corners of the frame.
  162.  */
  163.  
  164. void DrawFrame( int iTop, int iLeft, int iWidth, int iHeight ) {
  165.  
  166.     enum { ULEFT = 201, URIGHT = 187,
  167.            LLEFT = 200, LRIGHT = 188,
  168.            VERTICAL = 186, HORIZONTAL = 205
  169.            };
  170.     int iRow;
  171.     char achTmp[50];
  172.  
  173.     memset( achTmp, HORIZONTAL, iWidth );
  174.     achTmp[0] = (char)ULEFT;
  175.     achTmp[iWidth - 1] = (char)URIGHT;
  176.     achTmp[iWidth] = '\0';
  177.     _outxcr( achTmp, iLeft, iTop );
  178.  
  179.     memset( achTmp, BLANK, iWidth );
  180.     achTmp[0] = (char)VERTICAL;
  181.     achTmp[iWidth - 1] = (char)VERTICAL;
  182.     for(  iRow = iTop + 1; iRow <= iHeight; iRow++ )
  183.         _outxcr( achTmp, iLeft, iRow );
  184.  
  185.     memset( achTmp, HORIZONTAL, iWidth );
  186.     achTmp[0] = (char)LLEFT;
  187.     achTmp[iWidth - 1] = (char)LRIGHT;
  188.     _outxcr( achTmp, iLeft, iTop + iHeight );
  189.     }
  190.  
  191. int RunMenu() {
  192.  
  193.     char ch;
  194.  
  195.     if( dispMenu == 0 ) {
  196.  
  197.       if( !kbhit() )
  198.         return( 0 );
  199.  
  200.       else {
  201.         dispMenu = 1;
  202.         _clearscreen( _GCLEARSCREEN );
  203.         bwidth = 42;
  204.         Hbar[bwidth] = '\0';
  205.         fwidth = 42.5;
  206.         prevC = -1;
  207.         DrawMenu( MENUCOLOR );
  208.         sprintf( String, "<<<" );
  209.         if( dispType == 0 ) {
  210.             _outxcr( String, LEFTCOLUMN + 20, TOP + 4 );
  211.             }
  212.         else if( dispType == 2 ) {
  213.             _outxcr( String, LEFTCOLUMN + 20, TOP + 5 );
  214.             }
  215.         else {
  216.             _outxcr( String, LEFTCOLUMN + 20, TOP + 6 );
  217.             }
  218.         sprintf( String, "Replications =%10ld", reps );
  219.         _outxcr( String, LEFTCOLUMN, TOP + 12 );
  220.         sprintf( String, "Cum. Chi Sq. =%10.2lf (%2d)", chisqtot, doftot );
  221.         _outxcr( String, LEFTCOLUMN, TOP + 13 );
  222.         sprintf( String, "Inc. Chi Sq. =%10.2lf (%2d)", chisqinc, dofinc );
  223.         _outxcr( String, LEFTCOLUMN, TOP + 14 );
  224.         _displaycursor( _GCURSOROFF );
  225.         while( kbhit() )
  226.             getch();
  227.         return( 0 );
  228.         }
  229.       }
  230.  
  231.     _settextcolor( MENUCOLOR );
  232.     _setbkcolor( BACKCOLOR );
  233.     sprintf( String, "Replications =%10ld", reps );
  234.     _outxcr( String, LEFTCOLUMN, TOP + 12 );
  235.     sprintf( String, "Cum. Chi Sq. =%10.2lf (%2d)", chisqtot, doftot );
  236.     _outxcr( String, LEFTCOLUMN, TOP + 13 );
  237.     sprintf( String, "Inc. Chi Sq. =%10.2lf (%2d)", chisqinc, dofinc );
  238.     _outxcr( String, LEFTCOLUMN, TOP + 14 );
  239.  
  240.     if( !kbhit() )
  241.         return( 0 );
  242.  
  243.     _settextposition( TOP + cszMenu, LEFTCOLUMN + PROMPTPOS );
  244.     ch = (char)getch();
  245.     sprintf( MenuCh, "%c", toupper( ch ) );
  246.     _outtext( MenuCh );
  247.     while( kbhit() )
  248.         getch();
  249.     
  250.     sprintf( String, "   " );
  251.     _outxcr( String, LEFTCOLUMN + 20, TOP + 4 );
  252.     _outxcr( String, LEFTCOLUMN + 20, TOP + 5 );
  253.     _outxcr( String, LEFTCOLUMN + 20, TOP + 6 );
  254.     sprintf( String, "<<<" );
  255.  
  256.     /* Branch to the appropriate procedure depending on the key. */
  257.  
  258.     switch( toupper( ch ) ) {      /* choices = T C I L M P R S */
  259.         
  260.         case 'T':
  261.             dispType = 0;
  262.             break;
  263.     
  264.         case 'C':
  265.             dispType = 2;
  266.             break;
  267.     
  268.         case 'I':
  269.             dispType = 1;
  270.             break;
  271.     
  272.         case 'L':
  273.             if( tRows == 25 )
  274.                 tRows = 43;
  275.             else
  276.                 tRows = 25;
  277.             prevC = -1;
  278.             ClrScreen();
  279.             DrawMenu( MENUCOLOR );
  280.             break;
  281.  
  282.         case 'P':
  283.             _clearscreen( _GCLEARSCREEN );
  284.             SSSparm();
  285.             prevC = -1;
  286.             ResetHist();
  287.             ClrScreen();
  288.             DrawMenu( MENUCOLOR );
  289.             break;
  290.  
  291.         case 'R':
  292.             _clearscreen( _GCLEARSCREEN );
  293.             ResetHist();
  294.             prevC = -1;
  295.             ClrScreen();
  296.             DrawMenu( MENUCOLOR );
  297.             break;
  298.  
  299.         case 'S':
  300.             _clearscreen( _GCLEARSCREEN );
  301.             SampSize();
  302.             prevC = -1;
  303.             ClrScreen();
  304.             DrawMenu( MENUCOLOR );
  305.             break;
  306.  
  307.         case 'M':
  308.             dispMenu = 0;
  309.             prevC = -1;
  310.             Hbar[bwidth] = (char)BLOCK;
  311.             bwidth = 77;
  312.             fwidth = 77.5;
  313.             DrawMenu( MENUERASE );
  314.             break;
  315.  
  316.         case 'N':
  317.             if( cellwid == 2 )
  318.                 cellwid = 1;
  319.             if( cellwid == 4 )
  320.                 cellwid = 2;
  321.             prevC = -1;
  322.             break;
  323.  
  324.         case 'W':
  325.             if( cellwid == 2 )
  326.                 cellwid = 4;
  327.             if( cellwid == 1 )
  328.                 cellwid = 2;
  329.             prevC = -1;
  330.             break;
  331.  
  332.         case ESC:
  333.             return( 1 );
  334.     
  335.         default:        // Unknown key ignored
  336.             sprintf( MenuCh, " " );
  337.             break;
  338.         }
  339.  
  340.     _settextposition( TOP + cszMenu, LEFTCOLUMN + PROMPTPOS );
  341.     _outtext( MenuCh );
  342.  
  343.     if( dispType < 1 ) {
  344.             _outxcr( String, LEFTCOLUMN + 20, TOP + 4 );
  345.             }
  346.     else if( dispType > 1 ) {
  347.             _outxcr( String, LEFTCOLUMN + 20, TOP + 5 );
  348.             }
  349.     else {
  350.             _outxcr( String, LEFTCOLUMN + 20, TOP + 6 );
  351.             }
  352.  
  353.     _displaycursor( _GCURSOROFF );
  354.  
  355.     return( 0 );
  356.     }
  357.     
  358. void DrawHist( int Code ) {
  359.     int iRow, iLength, Separate, iColor;
  360.     double Fract;
  361.  
  362.     if( Code == 0 && prevC == 0 )
  363.         return;
  364.     prevC = Code;
  365.  
  366.     for( iRow = topRow; iRow <= botRow; iRow++ ) {
  367.  
  368.         Separate = 0;
  369.  
  370.         if( cellwid == 4 ) {
  371.             if( ( topRow ==  3 && iRow <  5 ) ||
  372.                 ( botRow == 45 && iRow > 44 ) )
  373.                 Separate = 1;
  374.             else if( ( topRow == 12 && iRow == 12 ) ||
  375.                      ( botRow == 36 && iRow == 36 ) )
  376.                 Separate = 1;
  377.             else if( iRow % 4 == 1 ) {
  378.                 if( Code == 0 )
  379.                     Fract = ( expfrac[iRow]   + expfrac[iRow+1] +
  380.                               expfrac[iRow+2] + expfrac[iRow+3] ) / 4;
  381.                 else if( Code == 1 )
  382.                     Fract = ( obsfrac[iRow]   + obsfrac[iRow+1] +
  383.                               obsfrac[iRow+2] + obsfrac[iRow+3] ) / 4;
  384.                 else
  385.                     Fract = ( cumfrac[iRow]   + cumfrac[iRow+1] +
  386.                               cumfrac[iRow+2] + cumfrac[iRow+3] ) / 4;
  387.                 iLength = (int)( EFRAC * fwidth * Fract / maxfrac );
  388.                 }
  389.             }
  390.  
  391.         if( cellwid == 2 ) {
  392.             if( ( topRow == 12 && iRow == 12 ) ||
  393.                 ( botRow == 36 && iRow == 36 ) )
  394.                 Separate = 1;
  395.             else if( iRow % 2 == 1 ) {
  396.                 if( Code == 0 )
  397.                     Fract = ( expfrac[iRow] + expfrac[iRow+1] ) / 2;
  398.                 else if( Code == 1 )
  399.                     Fract = ( obsfrac[iRow] + obsfrac[iRow+1] ) / 2;
  400.                 else
  401.                     Fract = ( cumfrac[iRow] + cumfrac[iRow+1] ) / 2;
  402.                 iLength = (int)( EFRAC * fwidth * Fract / maxfrac );
  403.                 }
  404.             }
  405.  
  406.         if( cellwid == 1 || Separate == 1 ) {
  407.             if( Code == 0 )
  408.                 iLength = (int)( EFRAC * fwidth * expfrac[iRow] / maxfrac );
  409.             else if( Code == 1 )
  410.                 iLength = (int)( EFRAC * fwidth * obsfrac[iRow] / maxfrac );
  411.             else
  412.                 iLength = (int)( EFRAC * fwidth * cumfrac[iRow] / maxfrac );
  413.             }
  414.  
  415.         if( iLength > bwidth )
  416.             iLength = bwidth;
  417.  
  418.         abar[iRow-topRow].length = iLength;
  419.         if( iColorMax == 1 )
  420.             abar[iRow-topRow].color = (short)SINGLCOLOR;
  421.         else {
  422.             iColor = (int)( (double)(iLength*14) / (double)bwidth + 1.0 );
  423.             abar[iRow-topRow].color = (short)iColor;
  424.             }
  425.  
  426.         DrawCell( iRow );
  427.         _displaycursor( _GCURSOROFF );
  428.         }
  429.     }
  430.  
  431. /*
  432.  * DrawCell - Prints a bar at a specified row using
  433.  * the length and color given in the abar work array.
  434.  */
  435.  
  436. void DrawCell( int iRow ) {
  437.     int    rValue, relRow;
  438.     char   Label[3];
  439.  
  440.     relRow = iRow-topRow;
  441.     _settextposition( relRow+1, 3 );        /* starting from column 3 */
  442.     _settextcolor( 0 );
  443.     _outtext( Hbar );                       /* erase 42 or 77 positions */
  444.  
  445.     if( abar[relRow].length > 0 ) {
  446.         _settextposition( relRow+1, 3 );    /* restart from column 3 */
  447.         _settextcolor( abar[relRow].color );
  448.         Hbar[abar[relRow].length] = '\0';
  449.         _outtext( Hbar );
  450.         if( abar[relRow].length < bwidth )
  451.             Hbar[abar[relRow].length] = (char)BLOCK;
  452.         }
  453.  
  454.     if( ( iRow ) % 4 == 0 ) {
  455.         if( dispMenu == 1 ) {
  456.             rValue = 6 - iRow / 4;
  457.             sprintf( Label, "%2d", rValue );
  458.             _settextcolor( SINGLCOLOR );
  459.             _outxcr( Label, 1, relRow+1 );
  460.             }
  461.         else {
  462.             memset( Label, BLANK, 2 );
  463.             Label[2] = '\0';
  464.             _outxcr( Label, 1, relRow+1 );
  465.             }
  466.         }
  467.     if( dispType == 1 && dispMenu == 1 && iRow <= firstobs ) {
  468.         sprintf( Label, "" );
  469.         _settextcolor( MENUCOLOR );
  470.         _outxcr( Label, 9, relRow+1 );
  471.         }
  472.     if( dispType == 1 && dispMenu == 1 && iRow >= lastobs ) {
  473.         sprintf( Label, "" );
  474.         _settextcolor( MENUCOLOR );
  475.         _outxcr( Label, 9, relRow+1 );
  476.         }
  477.     if( dispType == 2 && dispMenu == 1 && iRow <= firstcum ) {
  478.         sprintf( Label, "" );
  479.         _settextcolor( MENUCOLOR );
  480.         _outxcr( Label, 9, relRow+1 );
  481.         }
  482.     if( dispType == 2 && dispMenu == 1 && iRow >= lastcum ) {
  483.         sprintf( Label, "" );
  484.         _settextcolor( MENUCOLOR );
  485.         _outxcr( Label, 9, relRow+1 );
  486.         }
  487.  /* debugging: print bar lengths
  488.   * sprintf( Label, "%2d", abar[iRow].length );
  489.   * _settextcolor( SINGLCOLOR );
  490.   * _outxcr( Label, 12, iRow );
  491.   */
  492.     }
  493.  
  494. /*
  495.  * Standard Normal (Mean=0, Variance=1) Cumulative Distribution Function:
  496.  *
  497.  *      PHI( x ) = 1.0 - 0.5 * erfcc( x divided by square root of 2 )
  498.  */
  499.  
  500. void ClrScreen() {
  501.         int iRow, cRow;
  502.         cRow = _settextrows( tRows );   
  503.      /* cRow = number of screen rows actually available <= tRows */
  504.         if( cRow == 43 ) {
  505.                 topRow = 3;
  506.                 botRow = 45;
  507.                 }
  508.         else {
  509.              /* cRow   = 25; */
  510.                 topRow = 12;
  511.                 botRow = 36;
  512.                 }
  513.         _clearscreen( _GCLEARSCREEN );
  514.         _displaycursor( _GCURSOROFF );
  515.         rroot2 = 1.0 / sqrt( 2.0 );
  516.  
  517.         /* cutval = lower limit for a row and upper limit for next row down */
  518.         cutval[0]  = zvalue   = 6.0;
  519.         expfrac[0] = exptail[0] = pgreater = 0.5 * erfcc( zvalue * rroot2 );
  520.         for( iRow = 1; iRow < 49; iRow++ ) {
  521.                zvalue -= 0.25;
  522.                cutval[iRow] = zvalue;
  523.                pequal = 0.5 * erfcc( zvalue * rroot2 );
  524.                if( zvalue >= 0.0 )
  525.                    exptail[iRow] = pequal;
  526.                exptail[iRow+1] = 1.0 - pequal;
  527.                expfrac[iRow] = pequal - pgreater;
  528.                pgreater = pequal;
  529.                }
  530.         expfrac[49] = exptail[49];
  531.         cutval[49] = -99.99;
  532.  
  533.         maxfrac = 0.0;
  534.         for( iRow = topRow; iRow <= botRow; iRow++ )
  535.             if( expfrac[iRow] > maxfrac )
  536.                 maxfrac = expfrac[iRow];
  537.         }
  538.  
  539. void ResetHist() {
  540.       int iRow;
  541.  
  542.       if( iduminit == 0 )
  543.              idum = initseed();
  544.       else if( iduminit > 0 )
  545.              idum = -iduminit ;
  546.       else
  547.              idum =  iduminit ;
  548.  
  549.       if( savfile != NULL )
  550.              fprintf( savfile, "Random Number Seed  = %d\n", idum );
  551.  
  552.       reps = (long)0;
  553.       exptime( ltime );
  554.       for( iRow = 0; iRow < 50; iRow++ )
  555.            totals[iRow] = (long)0;
  556.       }
  557.  
  558. void ObsrvHist() {
  559.     int iRow, flagobs, flagcum;
  560.     double obstail, cumtail;
  561.  
  562.     for( iRow = 0; iRow < 50; iRow++ )
  563.          counts[iRow] = (long)0;
  564.  
  565.     for( iter = 0; iter < reppbat; iter++ ) {
  566.  
  567.         if( bmtype==0 )
  568.             zvalue = gasdev(&idum);
  569.         else
  570.             zvalue = box_mull(&idum);
  571.  
  572.         for( iRow = 0; iRow < 24; iRow++ ) {
  573.             if( zvalue < cutval[24-iRow-1] && zvalue >= cutval[24-iRow] ) {
  574.                 counts[24-iRow]++;
  575.                 break;
  576.                 }
  577.             if( zvalue < cutval[24+iRow] && zvalue >= cutval[24+iRow+1] ) {
  578.                 counts[24+iRow+1]++;
  579.                 break;
  580.                 }
  581.             }
  582.         if( zvalue >= 6.0 )
  583.             counts[0]++;
  584.         if( zvalue < -6.0 )
  585.             counts[49]++;
  586.         }
  587.     reps += (long)reppbat;
  588.     chisqtot = chisqinc = obstail = cumtail = 0.0;
  589.     dofinc = doftot = flagobs = flagcum  = 0;
  590.     firstobs = firstcum = -1;
  591.     lastobs = lastcum = 99;
  592.     for( iRow = 0; iRow < 50; iRow++ ) {
  593.          chisqobs[iRow] = chisqcum[iRow] = 0.0;
  594.          obsfrac[iRow] = (double)counts[iRow] / (double)reppbat;
  595.          totals[iRow] += counts[iRow];
  596.          cumfrac[iRow] = (double)totals[iRow] / (double)reps;
  597.  
  598.          if( flagobs == 1 && exptail[iRow+1] * (double)reppbat < exptmin ) {
  599.              flagobs = 0;
  600.              lastobs = iRow;
  601.              }
  602.          if( flagobs == 0 )
  603.              obstail += obsfrac[iRow];
  604.          else {
  605.              chisqobs[iRow] = square( obsfrac[iRow] - expfrac[iRow] ) / expfrac[iRow];
  606.              dofinc++;
  607.              }
  608.          if( flagobs == 0 && lastobs == 99 &&
  609.              exptail[iRow] * (double)reppbat >= exptmin ) {
  610.              chisqobs[iRow] = square( obstail - exptail[iRow] ) / exptail[iRow];
  611.              flagobs = 1;
  612.              firstobs = iRow;
  613.              obstail = 0.0;
  614.              }
  615.          chisqobs[iRow] *= (double)reppbat;
  616.          chisqinc += chisqobs[iRow];
  617.  
  618.          if( flagcum == 1 && exptail[iRow+1] * (double)reps < exptmin ) {
  619.              flagcum = 0;
  620.              lastcum = iRow;
  621.              }
  622.          if( flagcum == 0 )
  623.              cumtail += cumfrac[iRow];
  624.          else {
  625.              chisqcum[iRow] = square( cumfrac[iRow] - expfrac[iRow] ) / expfrac[iRow];
  626.              doftot++;
  627.              }
  628.          if( flagcum == 0 && lastcum == 99 &&
  629.              exptail[iRow] * (double)reps >= exptmin ) {
  630.              chisqcum[iRow] = square( cumtail - exptail[iRow] ) / exptail[iRow];
  631.              flagcum = 1;
  632.              firstcum = iRow;
  633.              cumtail = 0.0;
  634.              }
  635.          chisqcum[iRow] *= (double)reps;
  636.          chisqtot += chisqcum[iRow];
  637.          }
  638.  
  639.     if( lastobs > 0 ) {
  640.          chisqobs[lastobs] = square( obstail - exptail[lastobs] ) / exptail[lastobs];
  641.          chisqobs[lastobs] *= (double)reppbat;
  642.          chisqinc += chisqobs[iRow];
  643.          dofinc++;
  644.          }
  645.  
  646.     if( lastcum > 0 ) {
  647.          chisqcum[lastcum] = square( cumtail - exptail[lastcum] ) / exptail[lastcum];
  648.          chisqcum[lastcum] *= (double)reps;
  649.          chisqtot += chisqcum[lastobs];
  650.          doftot++;
  651.          }
  652.     }
  653.  
  654. /*
  655.  *      initseed()...initialize pseudo-random sequence using the Microsoft
  656.  *                   long time() function to read system clock.
  657.  */
  658.  
  659. int initseed() {
  660.     time_t stime;
  661.     int seed;
  662.  
  663.     time( &stime );
  664.     seed = (int)stime ;
  665.     if( seed > 0 ) return( -seed );
  666.     else if( seed == 0 ) return( -13 );
  667.     else return( seed );
  668.     }
  669.  
  670. /*
  671.  *      exptime()...time expired running SSS (using the Microsoft
  672.  *                   long time() function to read system clock.)
  673.  */
  674.  
  675. time_t exptime( ltime )
  676. time_t *ltime; {
  677.  
  678.     time_t etime;
  679.     static time_t begin = (long)0;
  680.  
  681.     time( &etime );
  682.     *ltime = etime ;
  683.     if( begin == (long)0 ) {
  684.             begin = etime ;
  685.             return (long)0;
  686.             }
  687.     else {
  688.             etime -= begin;
  689.             begin = *ltime;
  690.             return etime;
  691.             }
  692.     }
  693.  
  694. /*
  695.  *      ranX()...return rntype value from ran1(), ran2(), or ran3().
  696.  */
  697.  
  698. float ranX( idum )
  699. int *idum; {
  700.         extern rntype;
  701.         if( rntype == 1 )
  702.               return ran1( idum );
  703.         else if( rntype == 2 )
  704.               return ran2( idum );
  705.         else if( rntype == 3 )
  706.               return ran3( idum );
  707.         else
  708.               return whran( idum );
  709.         }
  710.  
  711. /*
  712.  *  Random number generator based upon: Wichman and Hill(1982),
  713.  *                                      Algorithm AS 183,
  714.  *                                      Applied Statistics 31, 188-190.
  715.  */
  716.  
  717. float whran( idum )
  718. int *idum; {
  719.         static int n1, n2, n3;
  720.         static double nn1, nn2, nn3;
  721.         static int ldum, ix=1, iy=1, iz=1; /* SEEDS */
  722.         static int iff=0;
  723.         double ixx, iyy, izz, rand;
  724.  
  725.         if( *idum < 0 || iff == 0 ) {
  726.                 ldum = (int)abs(*idum);
  727.                 *idum = 1;
  728.                 while( ldum > 30000 )
  729.                     ldum = ldum%(30000);
  730.                 if( ldum < 1 )
  731.                     ldum = 1;
  732.                 iff=1;
  733.                 nn1 = n1 = 30269;
  734.                 nn2 = n2 = 30307;
  735.                 nn3 = n3 = 30323;
  736.                 ix = iy = iz = ldum;
  737.                 }
  738.         do {
  739.             ix = 171 * (ix%177) -  2 * (ix/177);
  740.             if( ix < 0 )
  741.                 ix += 30269;
  742.             ixx = ix;
  743.  
  744.             iy = 172 * (iy%176) - 35 * (iy/176);
  745.             if( iy < 0 )
  746.                 iy += 30307;
  747.             iyy = iy;
  748.  
  749.             iz = 170 * (iz%178) - 63 * (iz/178);
  750.             if( iz < 0 )
  751.                 iz += 30323;
  752.             izz = iz;
  753.  
  754.             rand = fmod( ixx/nn1 + iyy/nn2 + izz/nn3, 1.0 );
  755.             } while( rand==0.0 || rand==1.0 );
  756.         return (float)rand;
  757.         }
  758.  
  759. void SSSparm() {
  760.  
  761.         printf( "\n\n\t***                   Six Sigma Simulator                   ***\n");
  762.         printf( "\n\n\t                    SSS.EXE....Version 9207" );
  763.         printf( "\n\n\t               A Quality Assurance Training Tool:" );
  764.         printf(   "\n\t       Statistics Committee of the QA Section of the PMA" );
  765.         printf( "\n\n\t           Bob Obenchain, CompuServe User [72007,467]\n\n" );
  766.  
  767.         strcpy( parnam, "sss.par" );
  768.  
  769.         printf(  "\n\n\tAt colon Prompts : ...simply press ENTER to get the [default]." );
  770.  
  771.         printf(  "\n\tThe filename for saving Keyboard input is to be : %s\n",
  772.                 parnam );
  773.  
  774.         printf(  "\n\n\tRespond to the following prompts for Parameters...\n" );
  775.  
  776.         printf(  "\n\n\tType of Box-Muller Transformation ?" );
  777.         printf(    "\n\tFast, Rejection Algorithm    (type=0) or" );
  778.         printf(    "\n\tSlow, Trignometric, Accurate (type=1)     [%d] : ",
  779.                 bmtype );
  780.         gets( buf );
  781.         if( buf[0] != '\0' && buf[0] != '\r' && buf[0] != '\n' )
  782.                 sscanf( buf, "%d", &bmtype );
  783.         if( bmtype > 0 )
  784.                 bmtype = 1 ;
  785.         else
  786.                 bmtype = 0 ;
  787.  
  788.         printf(  "\n\n\tType of Portable, Uniform Pseudo-Random Generator ?" );
  789.         printf(    "\n\t3 Linear Gongruential Gens. (type=1)" );
  790.         printf(    "\n\tFast, Single Generator COMB (type=2)" );
  791.         printf(    "\n\tSubtractive, Floating Knuth (type=3)  or" );
  792.         printf(    "\n\tWichman and Hill Algorithm  (type=4)     [%d] : ", rntype );
  793.         gets( buf );
  794.         if( buf[0] != '\0' && buf[0] != '\r' && buf[0] != '\n' )
  795.                 sscanf( buf, "%d", &rntype );
  796.         if( rntype > 3 )
  797.                 rntype = 4 ;
  798.         else if( rntype < 2 )
  799.                 rntype = 1 ;
  800.  
  801.         printf(  "\n\n\tStart-up Seed ( 0 => use Clock )   [%d] : ",
  802.                 iduminit );
  803.         gets( buf );
  804.         if( buf[0] != '\0' && buf[0] != '\r' && buf[0] != '\n' )
  805.                 sscanf( buf, "%d", &iduminit );
  806.         if( iduminit > 0 )
  807.                 iduminit = -iduminit ;
  808.  
  809.         printf(  "\n\n\tMinimum Expectancy (1 to 5) in Tail Cells [%4.1lf] : ", exptmin );
  810.         gets( buf );
  811.         if( buf[0] != '\0' && buf[0] != '\r' && buf[0] != '\n' )
  812.                 sscanf( buf, "%lf", &exptmin );
  813.         if( exptmin < 1.0 || exptmin > 5.0 ) {
  814.                 printf(  "\n\tMinimum Expectancy will be 3.0, the default.\n");
  815.                 exptmin = 3.0;
  816.                 }
  817.  
  818.         printf(  "\n\n\tMaximum Total Number of Replications [%ld] : ", maxreps );
  819.         gets( buf );
  820.         if( buf[0] != '\0' && buf[0] != '\r' && buf[0] != '\n' )
  821.                 sscanf( buf, "%ld", &maxreps );
  822.         if( maxreps > 2000000000 ) {
  823.                 printf(  "\n\tMaximum MAXREPS = 2,000,000,000.\n");
  824.                 maxreps = 2000000000;
  825.                 }
  826.         else if( maxreps < 10000 ) {
  827.                 printf(  "\n\tMinimum MAXREPS = 10,000.\n");
  828.                 maxreps = 10000;
  829.                 }
  830.  
  831.         if( ( parfile = fopen( parnam, "w" ) ) == NULL ) {
  832.                 printf(  "\tCannot write to Keyboard save file : %s\n",
  833.                         parnam );
  834.                 return;
  835.                 }
  836.  
  837.         fprintf( parfile, "%6d  ;number of replications per batch\n",
  838.             reppbat );
  839.         fprintf( parfile, "%6d  ;Box-Muller type [0,1]\n", bmtype );
  840.         fprintf( parfile, "%6d  ;random number type [1,2,3,4]\n", rntype );
  841.         fprintf( parfile, "%6d  ;random number seed\n", iduminit );
  842.         fprintf( parfile, "%6d  ;number of screen rows\n", tRows );
  843.         fprintf( parfile, "%6.1lf  ;minimum expectancy in tail cells\n", exptmin );
  844.         fprintf( parfile, "%6ld  ;maximum total replications\n", maxreps );
  845.         fclose( parfile );
  846.         }
  847.  
  848. void SampSize() {
  849.  
  850.         printf( "\n\n\t***                   Six Sigma Simulator                   ***\n");
  851.         printf( "\n\n\t                    SSS.EXE....Version 9207" );
  852.         printf( "\n\n\t               A Quality Assurance Training Tool:" );
  853.         printf(   "\n\t       Statistics Committee of the QA Section of the PMA" );
  854.         printf( "\n\n\t           Bob Obenchain, CompuServe User [72007,467]\n\n" );
  855.  
  856.         printf(  "\n\n\tSample Size now set at... %d Replications per Batch.",
  857.                 reppbat );
  858.  
  859.         printf(  "\n\n\tSpecify Sample Size (10 to 1000) [%d] : ", reppbat );
  860.         gets( buf );
  861.         if( buf[0] != '\0' && buf[0] != '\r' && buf[0] != '\n' )
  862.                 sscanf( buf, "%d", &reppbat );
  863.         if( reppbat < 10 || reppbat > 1000 ) {
  864.                 printf(  "\n\tNumber of Replications per Batch will be 100, the default.\n");
  865.                 reppbat = 100;
  866.                 }
  867.         }
  868.  
  869. void main() {
  870.  
  871.         int iRow;
  872.  
  873.         if( !_setvideomode( _TEXTC80 ) )
  874.             if( !_setvideomode( _TEXTMONO ) ) {
  875.                 _setvideomode( _DEFAULTMODE );
  876.                 printf( "\n\nYou cannot run SSS.  Your screen has only 40 Columns!" );
  877.                 exit( 0 );
  878.                 }
  879.  
  880.         /* If monochrome or color burst disabled, use one color */
  881.         _getvideoconfig( &vc );
  882.         if( (vc.monitor == _MONO) || (vc.mode == _TEXTBW80) )
  883.             iColorMax = 1;
  884.         else
  885.             iColorMax = 15;
  886.  
  887.         /* interrupt Ctrl-Break handler...return to text mode */
  888.         if( signal( SIGINT, handler ) == SIG_ERR ) {
  889.             printf("\nCouldn't set SIGINT...Abort!\n\n");
  890.             exit( 1 );
  891.             }
  892.  
  893.         strcpy( parnam, "sss.par" );
  894.         if( ( parfile = fopen( parnam, "r" ) ) == NULL ) {
  895.                 reppbat = 100;
  896.                 bmtype = 1;
  897.                 rntype = 1;
  898.                 iduminit = 0;
  899.                 tRows = 43; /* use 43 line mode if possible, else use 25 */
  900.                 exptmin = 3.0;
  901.                 maxreps = 2000000000; /* 2 billion */
  902.                 }
  903.         else {
  904.                 fgets( buf, 255, parfile );
  905.                 sscanf( buf, "%d", &reppbat );
  906.                 if( reppbat < 10 || reppbat > 1000 )
  907.                         reppbat = 100;
  908.                 fgets( buf, 255, parfile );
  909.                 sscanf( buf, "%d", &bmtype );
  910.                 if( bmtype > 0 )
  911.                         bmtype = 1 ;
  912.                 else
  913.                         bmtype = 0 ;
  914.                 fgets( buf, 255, parfile );
  915.                 sscanf( buf, "%d", &rntype );
  916.                 if( rntype > 3 )
  917.                         rntype = 4 ;
  918.                 else if( rntype < 2 )
  919.                         rntype = 1 ;
  920.                 fgets( buf, 255, parfile );
  921.                 sscanf( buf, "%d", &iduminit );
  922.                 fgets( buf, 255, parfile );
  923.                 sscanf( buf, "%d", &tRows );
  924.                 if( tRows != 25 )
  925.                         tRows = 43;
  926.                 fgets( buf, 255, parfile );
  927.                 sscanf( buf, "%lf", &exptmin );
  928.                 if( exptmin < 1.0 || exptmin > 5.0 )
  929.                         exptmin = 3.0;
  930.                 fgets( buf, 255, parfile );
  931.                 sscanf( buf, "%ld", &maxreps );
  932.                 if( maxreps > 2000000000 )
  933.                         maxreps = 2000000000;
  934.                 else if( maxreps < 10000 )
  935.                         maxreps = 10000;
  936.                 fclose( parfile );
  937.                 }
  938.  
  939.         if( iduminit == 0 )
  940.                 idum = initseed();
  941.         else if( iduminit > 0 )
  942.                 idum = -iduminit ;
  943.         else
  944.                 idum =  iduminit ;
  945.  
  946.         strcpy( savnam, "sss.sav" );
  947.         if( ( savfile = fopen( savnam, "w" ) ) == NULL )
  948.                 printf(  "\tCannot write to Save filename : %s\n", savnam );
  949.  
  950.         ClrScreen();
  951.         ResetHist();
  952.  
  953.         dispType = 1;   /* Incremental Histogram */
  954.         dispMenu = 0;   /* Menu Off */
  955.         prevC = dispType;
  956.         memset( Hbar, BLOCK, 77 );
  957.         bwidth = 77;
  958.         Hbar[bwidth] = '\0';
  959.         fwidth = 77.5;
  960.         cellwid = 1;
  961.  
  962.         while( TRUE ) {
  963.             ObsrvHist();
  964.             DrawHist( dispType );
  965.             if( RunMenu() == 1 || reps >= maxreps )
  966.                 break;
  967.             }
  968.  
  969.         if( savfile != NULL ) {
  970.             timer = exptime( ltime );
  971.             curtime = localtime( ltime );
  972.             fprintf( savfile, "Date/Time Stamp : %s", asctime( curtime ) );
  973.             fprintf( savfile, "Elapsed Simulation Time (sec) = %ld\n", timer );
  974.             fprintf( savfile, "Maximum Replications= %ld\n", maxreps );
  975.             fprintf( savfile, "\nSSS Parameter Settings...\n" );
  976.             fprintf( savfile, "%6d  ;Box-Muller type [0,1]\n", bmtype );
  977.             fprintf( savfile, "%6d  ;random number type [1,2,3,4]\n", rntype );
  978.             fprintf( savfile, "%6d  ;replications per batch\n",
  979.                 reppbat );
  980.             fprintf( savfile, "%6d  ;number of screen rows\n", tRows );
  981.             fprintf( savfile, "\nSSS Summary Statistics...\n" );
  982.             fprintf( savfile, "Total Number of Replications = %ld\n", reps );
  983.             fprintf( savfile, "Minimum Expectancy in Tail Cells = %4.1lf\n", exptmin );
  984.             fprintf( savfile, "Cum. Chi Sq. = %8.2lf with d.f.= %2d\n", chisqtot, doftot );
  985.             fprintf( savfile, "Inc. Chi Sq. = %8.2lf with d.f.= %2d\n", chisqinc, dofinc );
  986.             fprintf( savfile, "\nRow  CutVal    ExpFract     ExpTail    CumFract   ChiSq    ObsFract   ChiSq\n" );
  987.             for( iRow = 0; iRow < 50; iRow++ )
  988.                 fprintf( savfile, "%2d %8.4lf %11.9lf %11.9lf %11.9lf %7.2lf %11.9lf %7.2lf\n",
  989.                     iRow, cutval[iRow], expfrac[iRow], exptail[iRow],
  990.                     cumfrac[iRow], chisqcum[iRow], obsfrac[iRow], chisqobs[iRow] );
  991.             }
  992.  
  993.         _setvideomode( _DEFAULTMODE );
  994.         exit( 0 );
  995.         }
  996.  
  997.